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Abstract. Why life persists at the edge of chaos is a question at the very heart of evolution. 
Here we show that molecules taking part in biochemical processes from small molecules to 
proteins are critical quantum mechanically. Electronic Hamiltonians of biomolecules are tuned 
exactly to the critical point of the metal-insulator transition separating the Anderson localized 
insulator phase from the conducting disordered metal phase. Using tools from Random 
Matrix Theory we confirm that the energy level statistics of these biomolecules show the 
universal transitional distribution of the metal-insulator critical point and the wave functions 
are multifractals in accordance with the theory of Anderson transitions. The findings point 
to the existence of a universal mechanism of charge transport in living matter. The revealed 
bio-conductor material is neither a metal nor an insulator but a new quantum critical material 
which can exist only in highly evolved systems and has unique material properties. 


Advances in the theory of complex systems over the last quarter century reinforced that 
living systems exist at the edge of chaos [H [2f [3, 5] and order, poised at criticality[5j. Finding 
the detailed mechanism behind this apparent self-organized criticality pal .7] is still a tantalizing 
problem. One of the fascinating aspects of life is the highly organized molecular machinery taking 
care of myriads of complex processes such as DNA replication, protein synthesis, cell division 
and metabolism, to mention only a few. Electric forces animating the parts require a perpetual 
and precise motion of charges throughout the system for perfect execution of biochemical tasks. 
In this paper we show that practically all bionrolecules, from small signalling molecules to 
proteins taking part in biochemical electronic processes, belong to a fundamentally new class of 
conducting material. This is a disordered conductor where the strength of the disorder is tuned 
exactly to the metal-insulator transition point and it is consequently in a permanent critical 
quantum state. 

Our initial perspective is that of condensed-matter physics. The unique properties of the 
critical quantum state at the localization-delocalization transition point have been described 
first in the the well known Anderson[8] model, which represents the current paradigm for 
understanding conduction in condensed matter. In the Anderson Hamiltonian H = J2j e i a t a * — 
a+a; on a 3D lattice with random uniformly distributed on-site energies g E [—IT, IT], 
the critical level of disorder[5] is W c = 16 ± 0.5. For IT < IT C the system is a disordered 


metal with extended states and for W > W c the states are localized and the system is an 
insulator. It has been shown [TU] that in the exact critical point of the Anderson transition (W = 
W c ) the electron eigenfunctions are extended, but strongly inhomogeneous multifractals[11]. 
Similar metal-insulator transitions (MIT) exist in a wide range of physical systems at various 
dimensionalitiespL2lj including the ID Harper-Hofstadter model[13j and quantum-Hall-type 
transitions. 

For example, the ID Harper modelp3j describes the energy spectrum of an electron in 
a 2D lattice, whose graphic representation is widely known as the Hofstadter butterfly)13j. 
In the related Aubry-Andre HamiltonianjT5] H = ^ cos(2ircri)af a,i — where 

cr = (\/5 — l)/2) is the Golden Mean Ratio, the critical value separating localized and extended 
phases is A c = 2. The wave functions are also multifractals(l6j there. Other examples includef!2l 
quantum-Hall-type transitions in disordered conductors and superconductors in strong magnetic 
fields, the spin quantum Hall effect, the thermal quantum Hall effect, Dirac fermions in random 
vector potentials and Bethe lattices. 

MIT-like transitions exist also in low-dimensional quantum chaos. Quantum counterparts 
of strongly chaotic systems[TT] share the properties of delocalized systems, while integrable 
systems [18] have localized wave functions in quantum numbers corresponding to conserved 
quantities. Pseudointegrable systems[19j he at the border of chaos and integrability, where 
classical trajectories diverge only slowly (with zero Lyapunov exponent) but their dynamics 
is complex and the periodic orbits proliferate in their phasespace exponentially. Quantized 
pseudointegrable systems also show all the key features of critical systems, including multifractal 

wavefunctions[2Dj • 

Criticality can also be observed in the energy spectrum of systems at the transition point. 
Random Matrix Theory [21, 1221 23l [24] (RMT) is the main tool for the characterization of 
the universal statistical properties of Hamiltonians of complex systems. The distance between 
consecutive energy levels fluctuates in the spectrum. The raw distance between levels a n = 
E n+ i — E n can be normalized using the average separation of levels A (E) at a given energy 
window around E. The ratio s n = a n /A(E n ) is called the level spacing. Random matrix theory 
has certain predictions for the form of the distribution P(s) of level spacings. 

It has been established[25] that in the 3D Anderson model the localized, delocalized and 
the critical states each have a distinct level spacing distribution. These three distributions are 
believed to be universal, i.e. independent of the microscopic details of the disordered system. 
In the delocalized metallic phase (W < W c ) the distribution coincides with the level statistics 
of the Gaussian Orthogonal Ensemble (GOE), which is the ensemble of real symmetric random 
matrices with identically distributed Gaussian elements. The level spacing distribution is the 
Wigner surmise [21] 

Pw{s) = y exp 

In the localized insulating phase {W > W c ) the energy levels form a random Poisson process 
and the level spacing distribution is exponential 

Pp(s ) = exp (-s). 

It has been shown (T7, .12, 26| that not only in the 3D Anderson model, but in all other examples 
of Anderson-like transitions, the energy level statistics in the delocalized phase is universal and 
corresponds to the proper Gaussian Random Matrix Ensemble reflecting the symmetries of the 
system, while in the localized phase it is always random Poissonian|T8[ 12, 26]. 

At the transitional point {W = W c ) a third kind of intermediate spectral statistics Pt(s) 
exists, which is the hallmark of the critical state[25j|. In solid state models theoretical 



arguments [571[28] and numerical studies[25;, 3Dj suggest the general form 

Pt(s) = cisexp (-c 2 s 1+7 ) , 

where ci and c 2 are 7 dependent normalization constants. In the case of the 3D Anderson 
mo del [5U1| 7 ~ 0.2. In the Harper model [ 25 ] and in other pseudointegrable models [3Tj showing 
critical quantum chaos a numerical value of 7 ~ 0 has been found [ 25 ]. which supports the 
semi-Poissonian distribution 


P T (s) = Psp(s) = 4s exp (-2s). 

The semi-Poissonian distribution has also been predicted from a short range plasma model [321 
of energy levels introduced in RMT. 

While in physical systems the critical state can be reached only upon a careful tuning of the 
strength of the disorder, in the following we show that certain biomolecules are precisely at the 
critical state without any external tuning. 

The exact numerical solution of the Schrodinger equation for the electronic states of large 
molecules such as proteins is a prohibitive task. Various approximations have been developed, 
which reduce the problem to a one-electron problem in the effective field of the remaining 
electrons. Wave functions of molecules are usually written in the form of Linear Combinations 
of Atomic Orbitals (LCAO) fa = J2 r CirXr where fa is a Molecular Orbital (MO) represented 
as the sum of atomic orbital (AO) contributions Xr- The one-electron eigenenergies and 
eigenvectors then can be determined from the generalized eigenvalue equation HC = ESC , 
where S rs = (x r \ Xs) and H rs = (y> | H e ff \ Xs ) are the overlap and effective Hamiltonian 
matrices respectively. The effective Hamiltonian depends on the coefficients which makes the 
problem nonlinear in C. This is the case in Hartree-Fock and Density Functional Theory (DFT) 
calculations, which then cannot be routinely carried out for proteins involving thousands of 
atoms. If we restrict our interest to the localization-delocalization problem in valence electrons 
and treat the two-electron part of the Hamiltonian as in the case of the electrons in metals, 
in an average sense only, we can apply semi-empirical methods. Once the positions of the 
atoms are known the Extended Hiickel (EH) Molecular Orbital Method[33] is quite successful 
in calculating the MOs of organic molecules. The diagonal part Hfa ; of the EH Hamiltonian 
is given by the ionization energies of the AOs [34]. while the off-diagonal elements are calculated 
from the diagonal elements and the overlap matrix 

Hji f H) = -K ( H rr + H ss ) S rs , 

where the common choice for the empirical constant is K = 1.75. This is similar in spirit to 
other tight binding Hamiltonians in various models of the Anderson transition. 

Myoglobin is the first[35] and one of the best studied protein structures. It plays a central 
role in the oxygen storage of muscles. It consists of 153 amino acids and weighs about 18000 
Daltons. For the numerical studies we selected NMR data of a solution form|36] (PDB ID:1MYF) 
from the Protein Data Bank [37] of RCSB, as it is captured in the ’’living” state and contains 
the positions of hydrogen atoms essential for the calculations. The EH calculations have been 
carried out by the numerical package YAeHMOP. There are N = 6329 valence electron AOs 
and the EH Hamiltonian and overlap matrices are sufficiently large of dimension 6329 x 6329, 
which makes it possible to make a good numerical comparison with similar calculations in 
solid state physics[16] 113011251 23] and quantum chaos[26, 20. [3T| 123 ] , Lowdin transformation of 
the coefficient vector C = S l ^ 2 C has been applied to transform the EH Hamiltonian to a real 
symmetric self-adjoint operator H L = S' -1 / 2 H^ eh ^ S ' -1 / 2 , which satisfies a normal eigenequation 
H l C l = EC l . 


The Highest Occupied Molecular Orbital (HOMO) and the Lowest Unoccupied Molecular 
Orbital (LUMO) play a key role both in electron transport and reactions of organic molecules.For 
the visual demonstration of the fractal nature of the eigenfunctions we show the HOMO and the 
LUMO of Myoglobin in Fig. 1. Absolute values of coefficients are shown such that the index 
r is ordered in the sequence of appearance of atoms and orbitals in the amino acid sequence 
of the protein. Statistical similarity of the magnified part of the wave functions to the entire 
function is a visual indication of a fractal. 

The multifractal analysis of the protein wave functions is based on the standard box counting 
procedure|38|. dividing the ID index space of ordered AO indices along the protein sequence 
into Ni ~ N/l boxes of size l and determining the box probability of the wave function in the 
kth. box, 

i -1 

= E REW ^)| 2 ,k = l,...,N h 

n =o 

as a suitable measaure. If the gth moments of this measure are counted in all boxes and is 
proportional to some power r(q) of the box size, 

d) 

\ k / E 

multifractal behaviour might be derived. For a simple monofractal r(q) = (q — 1 )D, where D is 
the fractal dimension. For multifractals the r{q) curve is nonlinear and the generalized fractal 
dimensions D q can be recovered r(q) = (q — 1 )D q = lim^o Fix?/InL 

We expect that for extended wave functions in the conducting phase the coefficients are 
evenly distributed around their mean, which is C^\ 2 ^ ^ = (1/IV) Ylk=i \Cr(Ek)\ 2 = 1/iV, and 
are independent of the position due to the normalization of the wave functions. The measure 
then scales like Hk{l) ~ l/N and the moments scale as \q ( l/N)( q yielding D q = D = 1 
independent of q. For localized states in the insulating phase the coefficients are nearly zero 
except in a short interval of the size of the localization length £. The localization length is 
much smaller than the system size (<JV, therefore at intermediate length scales ^ < JV an 

interval of size l either contains the localization interval and almost the full probability fj,k(l) ~ 1 
or it is practically empty Hk(l) ~ 0. The moments do not scale with the length Xq ~ 1 and 
r(q) = 0 yielding D q = D = 0. 

The numerical values of the fractal dimension of the protein wave function can be determined 
more conveniently by calculating the box probability for all possible boxes of length l and 
performing an additional averaging to smooth out statistical fluctuations 

(? - N ‘ X ^ W,)) ° "Tirr, 1' ( (S ' c -' 2 ) 1 ■ 

In Fig. 2. we show the generalized fractal dimensions for Myoglobin obtained numerically. 
The most significant value is the correlation dimension D 2 ~ 0.5 which is just midway between 
localization D = 0 and delocalization D = 1 confirming that the system is critical and the wave 
functions are multifractals. We note that the same numerical value D 2 = 0.5 has been obtained 
also for critical quantum chaos[26, 20 ] , 

Next we show that the level statistics of Myoglobin is also transitional. The energy levels 
computed with the EH method have been analyzed with a statistical method of RMT suitable 
for the analysis of a relatively low number of eigenvalues. The distance between two consecutive 
levels < 7 i = Ei+i — Ei is normalized with the average of k level spacings to the left and to the 
right Aj = 2 FPI giving the unfolded level spacing s* = a l /A l = (2 k + 1 )(£) +1 - 





Figure 1. The HOMO/LUMO orbitals for Myoglobin (PDB ID:1MYF) calculated 
with the Extended Hiickel method. Vertical axis: the absolute value of C ] r in red for the 
HOMO and LUMO in green (flipped). Horizontal axis: the index sequence r of AOs ordered 
along the amino acid sequence. Inset: enlarged part of the box. 


Ei)/(Ei_ |_fc + i — Ei_k). The choice of k depends on the variability of the density of energy levels. 
Statistical averaging would require large k values, while fast variation (especially singularities) 
in the density of energy levels restricts our choice to low values. We have found that a choice 
of k = 2 ... 5 ensures the stability of the distribution in our examples. Next, following standard 
procedures[31]. the cummulative spacing I(S) = #{sj < S}/N = / 0 S P(s)ds is calculated and 
compared to the theoretical predictions. For the Poissonian statistics Ip(S ) = 1 — exp (—S'), 
for the Wigner surmise Iw{S) = 1 — exp (—7rS 2 /4) and for the semi-Poissonian transitional 
statistics Isp{S) = 1 — (2 S +1) exp (—2S). In the main part of Fig. 3. we show the cummulative 
level spacing for the 6328 spacings in the spectrum of Myoglobin and in the inset we show the 
difference to Isp(S). We can see that without any parameter fitting the spacings for Myoglobin 
follow the critical theoretical curve with astonishing precision. Note, that no parameter fitting 
is involved in the procedure, the calculated spacing distribution has a less than 3% error like 
in the case of systems of critical quantum chaos [26i [ 20 ]. which are purely theoretical models as 
opposed to our case, where the positions of atoms in the protein are obtained experimentally. 

In addition to Myoglobin, we included in Fig. 3. two other proteins with known biochemical 
functions, selected randomly from the PDB just by size (close to 10000 valence AOs) and 
by the availability of the necessary structural data including the coordinates of Hydrogen 
atoms. Human profilin (PDB ID:1PFL)[39| in solution form has N = 5232 valence AOs. 
It is a ubiquitous eukaryotic protein that binds to both cytosolic actin and the phospholipid 
phosphatidylinositol-4,5-bisphosphate. Human apolipoprotein E (PDB ID:2L7B)|l0] has N = 
11980 valence AOs. It is one of the major determinants in lipid transport, playing a critical role 
in atherosclerosis and other diseases. One can see that these other randomly picked proteins are 
also on the critical curve with the same precision as Myoglobin. Further analysis (not shown 
here) reveals essentially the same D 2 ~ 0.5 values and the same generalized dimension D q 
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Figure 2. Generalized fractal dimensions D q of wave functions of the protein 
Myoglobin (PDB ID:1MYF) averaged for all energies from the Extended Hiickel 
calculation. Inset: Scaling of \q f° r Q = 2 as a function of the box length l on a double 
logarithmic plot. The green linear function is fitted to the scaling region and its slope yields the 
correlation dimension D 2 = 0.5 ± 0.01. 
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spectra for these proteins as well. 

We should emphasize again, that finding a large tight binding Hamiltonian tuned exactly or 
almost exactly to the critical point by random chance can happen only with an astronomically 
low probability. So, finding just a single protein with more than 100 amino acids having this 
property at random is impossible. 

Next, we investigate whether criticality is restricted to certain proteins only or is it a more 
wide spread phenomenon. The verification of fractality is not possible for smaller molecules as 
it requires a length scale of two decades (say l ~ 10 — 1000) to fit a reliable power law to the 
curve x q (0 rs j l T ( q \ In the case of level statistics we need much less data to verify the shape 
of the cumulative level spacing 7(5). In the RMT analysis in solid state physics and quantum 
chaos normally on the order of a thousand levels is sufficient. Here we have found that using 
the technique developed for smaller data sets and the usage of cumulative level statistics 7(5) 
instead of the distribution function P(s), which requires the binning of the data, jointly allow us 
to verify molecules with as low as N = 80 valence AOs reliably. In Fig. 4. we show a collection 
of level statistics coming from organic molecules of various size. The 3D structures of small 
molecules were taken from PubChem|41| and the energies are calculated with the EH method 
(YAeHMOP). 

Our first and most striking observation is that each molecule investigated from the biological 
domain belongs to one of the ’’clean” categories Ip,I\y or Isp- The reason is by no means 
obvious. In RMT|23:j and especially in quantum chaos[H2j intermediate distributions between 
the classes occur in finite systems. It seems that the Hamiltonians of these biomolecules are not 
random, they are tuned firmly to one of the classes. 

Our first example is Silk (PDB ID:1SLK)[33], a protein which serves as a structural material 








Figure 3. Cummulative level spacing distribution I(S) for Myoglobin (blue triangles), 
Profilin (red circles) and Human Apolipoprotein E (green circles) are shown. Only every 300 
value is shown for legibility. For comparison the theoretical curves Ip(S ) (red line), Iw(S) (blue 
line) and Isp(S) (black line) are also shown. Inset: Difference between the data curves and 
Isp(S) are shown with the colour of the data set. Each value is plotted. The error is below 3% 
in probability. 


and does not take an active part in biochemical processes. In Fig. 4. we can see that it has 
Poissonian level statistics and it belongs to the localized class. This is in line with the fact 
that silk is a very good insulator. This example confirms that criticality is not the property 
of individual amino acids. Amino acid sequences can produce not just critical materials but 
insulators as well. 

Other examples for structural biomaterials in Fig. 4. are Dextrin (CID 62698), which is 
a gum like substance and Octadecane (CID 11635) which is an alkane hydrocarbon found in 
mineral oil and Gefarnate (CID 5282182) a water insoluble terpene fatty acid. They all show 
Poissonian level statistics and are good insulators. 

It is less obvious why DNA belongs to this category. In Fig. 4. we show level statistics for a 
21 basis pair DNA sequence (NDB ID:2JYK)filj which is clearly Poissonian. 

The conductivity properties of DNA or RNA sequences are highly debated, however there 
seems to be a consensus that native DNA is a wide band gap semiconductor, practically an 
insulator [35]. 

We picked a few representatives from essential classes of biomolecules. From each class we 
used the molecule with the largest number of AOs possible among all possible molecules having 
complete 3D crystallographic data in PubChem. They all show criticality and semi-Poissonian 
statistics. The list includes Linoleic acid (CID 5280450), Primary fluorescent chlorophyll 
catabolite (CID 54740347), Sucrose (CID 5988), Vitamin D3 (CID 25245915), Vitamin B12 
(CID 16212801) and the largest amino acid Leucine (CID 6106). We investigated a dozen 
more biomolecules of lesser size. For all practical purposes they all showed semi-Poissonian 
statistics, but the statistics had larger errors due to the small number of levels. This list includes 



















Figure 4. Cummulative level spacing distribution for various molecules. (Only 10-15 
data points for each molecule is shown for legibility.) The main part: Molecules with critical 
cummulative level statistics Isp(S) = 1 — (2 S + l)e“ 2S ) The molecules shown are Vitamin B12, 
Vitamin D3, Linoleic Acid, Primary fluorescent chlorophyll catabolite, Sucrose and Leucine. 
Middle part: Molecules with Poissonian cummulative level statistics Ip(S) = 1 — e~ s (Notice, 
that the inset is zoomed and the main axis scales don’t apply). The molecules shown are Dextrin, 
Silk, Octadecane, Gefarnate and a 21 base pair DNA sequence (NDB ID: 2JYK). Small part: 
Molecules with Wigner cummulative level statistics I\v{S) = 1 — e -71 " 5 "/ 4 (Notice, that the 
inset is zoomed and the main axis scales don’t apply). The molecules shown are Testosterone, 
Progesterone, Dibenzo(a,e)pyrene and Aristolochic Acid. 


Nicotine, Adenosine, Caffeine, Amphetamine, Benzoanthracene, Chlorpozamine, Glucose, Fatty 
acids omega 3 and omega 6, Picrotin, Picrotoxin, Theophylline, Thiactin, Xanthine. We also 
carried out the analysis for the 20 amino acids coded by the universal genetic code. They are 
generally too small for the analysis of the level statistics to the level shown in Fig. 4. but some 
of them can already be classified by inspection. Based on this it is likely that Arginine, Cysteine, 
Sele Histidine, Isoleucine, Leucine, Methionine, Phenylalanine, Proline, Serine, Threonine 
and Tryptophan are critical and Alanine, Asparagine, Aspartic-acid, Glutamic-acid, Lysine, 
Tyrosine and Valine show Poissonian level statistics. For Glutamine and Glycine the results are 
inconclusive. 

We can also find molecules which belong to the good conductor class with Wigner level 
spacing statistics. We could find this only in polycyclic molecules with delocalized wave 
functions spreading through the entire molecule. In Fig. 4. we show Testosterone (CID 6013), 
Progesterone (CID 5994), Dibenzo(a,e)pyrene (CID 9126) and Aristolochicacid (CID 2236). 
Among bionrolecules we could find only Steroids in this class, the rest of such molecules were 
involved in combustion such as polycyclic aromatic hydrocarbons and were toxic or carcinogenic. 

We can summarize these findings as follows: Most of the molecules taking part actively 
in biochemical processes are tuned exactly to the transition point and are critical conductors. 






There is only the special class of polycyclic molecules with closely packed aromatic rings which 
show metallic behaviour and delocalization and the class of ’’structural materials” which play a 
role in the mechanical stiffness of biological systems. Individual amino acids can be Poissonian 
or critical but they seem to form proteins and polypeptides which are also either Poissonian or 
critical. 

These findings suggest an entirely new and universal mechanism of conductance in biology 
very different from the one used in electrical circuits. In metallic conductors charges float due 
to voltage differences. The electrical field accelerates electrons while scattering on impurities 
dissipates their energy fixing a constant average propagation velocity. In biological systems 
we seldom see examples for this. A more likely scenario is that a charge entering a critical 
conductor biomolecule will be under the joint influence of the quantum Hamiltonian and the 
excessive decoherence caused by the environment|46]. Such conductance mechanism has been 
found for the excitons in light harvesting systems |47] and it is currently in the focus of research in 
Quantum Biology [381. In these systems Environment-assisted Quantum Transport [T9] (ENAQT) 
is dominant and facilitates the fast quantum spreading of excitations over the system. We think 
that this mechanism is more universal in biological systems and charges in biological conductors 
are also subjects of this transport mechanism. Recently we have shown[50] that ENAQT is the 
most effective at the critical point of the localization-delocalization transition and the excitonic 
Hamiltonians of light harvesting systems are also at or near the critical point |51]. In the localized 
regime transport is hindered by strong quantum effects, while in delocalized systems decoherence 
destroys quantum propagation and the anti-Zeno effect[19} slows down diffusion. At the mobility 
edge the existence of extended multifractal wave functions throughout the system ensure end- 
to-end transport while coherence decays only algebraically |50j ensuring a longer coherence time 
and supressing the anti-Zeno effect. 

Our results also suggest that quantum transport played a distinguished role in evolution and 
selection. Both the number of known small molecules and proteins is about 10 8 and the number 
of chemically feasible small (< 500-Da) organic compounds is astronomical, estimated|52j to 
be 10 60 . The number of proteins grows exponentially with the number n of amino acids as 
~ 20 n , and the largest known has about n ~ 26000. This shows that chemical and biological 
evolution selected only a tiny fraction p ~ 10~ 50 of possible small biomolecules and even less 
for proteins. As the probability of finding a critical molecule or protein by random chance 
is also astronomically low, the large number of critical molecules and proteins found by quasi 
random browsing of major databases just by size and availability of the 3D crystallographic data 
suggests that criticality of the quantum Hamiltonian is prevailing in the evolutionary selection 
of biomolecules. Besides, the fact that some proteins are natural critical conductors may open 
up new avenues in materials science as well. 
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